!p.multi=[0,2,1]
!p.charsize=1.5
l=64
rr=1.      ;m
rmin=0.
rmax=2.
rds=rr     ;m
rdsi=1./rds
dz=(rmax-rmin)*rr/double(l-1)
g=10.    ;m/s^2
g0=g
z = indgen(l)*dz

; bottom values

rho00=1. 
p00=8
T00=20
th00=22. 

cp=1.
Rg=cp/2.5
cap=rg/cp
capi=1./cap
ss=g0/(cp*th00)

;rho = rho00*(1.-(ss*z)/(1.+z*rdsi))^(capi-1)
rho=rho00*(1.-ss*z)^(capi-1.)

plot,z,rho

st=-3./(th00*(l-1)*dz)
the=th00*(1.+ st*z)
plot,r/rr,the,yr=[min(the),max(the)]

!p.multi=0
end

